clc
clear all
L =10;
N  = 128;
n2 = N/2;
S = peaks(N);% rand(128,128);% full(createPsiK2D([8 0],[0.05 0.05],[0 0],L,N));% [1:128]'* [1:128];% peaks(N);
filterLen = 10;
N2 = N/2;


for k = -100:100
[Fm W W2] = creatingFilterMatrix(N,filterLen,k);
%fourier wavelet filter
Fnh = createFn(N2);
blkdiag(Fnh);
fnd = blkdiag(Fnh,Fnh);
A = diag(W2(1:N2,1));
C = diag(W2((N2+1):end,1));
B = diag(W2(1:N2,end));
D = diag(W2((N2+1):end,end));

A2 = (W2(1:N2,1));
C2 = (W2((N2+1):end,1));
B2 = (W2(1:N2,end));
D2 = (W2((N2+1):end,end));
Fw = [A B;C D];
%Wavelet transformation 

ac = (W * S.').';
ws = W * ac;
if(norm(ws-wt)< 1e-8)
    k
end
subplot(3,1,1)
imagesc(wt)
subplot(3,1,2)
imagesc(ws)

subplot(3,1,3)
imagesc(ws-wt)
drawnow
end

%fourier
t = fft2(S);
tic
t3 = fft(fft(S.').');
toc
tic
tb = myFFT(myFFT(S.').');
toc
tic
ac = (W * S.').';
ws = W * ac;
fc = ((Fw * fnd) * ws.').';
fs = (Fw * fnd) * fc;
toc

fo = (Fm * S.').';
fos = Fm * fo;
subplot(3,1,1)
imagesc(abs(fos))
subplot(3,1,2)
imagesc(abs(t3))
title('fft2')
subplot(3,1,3)
imagesc(abs(fs))
title('faft')
max(fs(:))


mwt = Fw * fnd*(Fw * fnd * ws.').';
mwt2 = Fw *(Fw * ( fnd * (fnd * ws.').').').';
norm(t3-fs)

w =  dbaux(filterLen);
[lpd,hpd,lpr,hpr] = orthfilt(w); 
x = (1:128)';
a = fft(x);
b = W*x;
[lp hp] = dwt(x,lpd,hpd,'mode','per');
b = [lp ; hp];
c = [fft(b(1:n2)); fft(b((n2+1):end))] ;
c2 = fnd*b;
subplot(3,1,2);
plot(imag(a))

%for k = -100:100
    [ A2 B2 C2 D2 ] = creatingCircRecon( N,filterLen,-filterLen +1 );
    d =  [(A2.*c(1:n2) + B2.*c((n2+1):end)) ;(C2.*c(1:n2) + D2.*c((n2+1):end))];
    subplot(3,1,3);
    plot(imag(d))
    subplot(3,1,1);
    plot(imag(d)-imag(a))
 

    
%end
    t = Fw*c;

norm(a/norm(a) - d/norm(d))

subplot(3,1,2);
plot(real(a))
subplot(3,1,3);
plot(real(d))
subplot(3,1,1);
plot(real(d)-real(a))


subplot(3,1,2);
plot(imag(a))
subplot(3,1,3);
plot(imag(d))
subplot(3,1,1);
plot(imag(d)-imag(a))

%subplot(3,1,1)
%plot(b)

